Confirm the best fit of a stochastic fit by running the top parameter sets again (#659) - #675

Closes #659.

The problem

When a model is stochastic, running it twice with the same parameter values gives two different answers. The objective value PyBNF computes is therefore a noisy measurement rather than a fixed number.

A fit picked its answer by taking the best objective value it ever saw, and every one of those values came from a single simulation. A long fit scores tens of thousands of parameter sets, so the winner of that comparison was very often just the parameter set that happened to get a lucky simulation.

Two things came out wrong and nothing in the output said so. The reported objective value was the best of many noisy draws, so it was optimistic by a wide margin. And the reported parameter values were not the best ones found, because a slightly worse parameter set with a lucky draw beats a better one with an average draw.

The change

A fit that uses at least one stochastic model now ends with one more stage. It takes the top best_fit_candidates parameter sets, runs each of them best_fit_replicates more times, ranks them by their average objective value, and reports that winner as the best fit.

Ten and ten by default. That is a hundred simulations against a run that may have done a hundred thousand, and all of them are submitted at once, which uses processors that would otherwise sit idle at the end of a run.

Results/best_fit_confirmation.txt records the result. Here is the real file from the run described below:

candidates	3
replicates_requested	4
winner	iter2p0h4
winner_mean_objective	72838.74908
winner_standard_error	5521.750431
winner_search_objective	54877.73818
optimism	17961.0109
# rank	name	mean_objective	standard_error	std_deviation	runs	failed	search_objective
1	iter2p0h4	72838.74908	5521.750431	11043.50086	4	0	54877.73818
2	init4	87262.13928	3134.119241	6268.238482	4	0	99030.06578
3	iter2p5h3	100953.6321	3221.85557	6443.71114	4	0	101437.2518

The winner is recorded on the trajectory rather than passed around, so everything the run writes after this stage describes the same parameter set. That includes the saved simulations, the best-fit model file, the information criteria, a refine's starting point, and a bootstrap replicate's recorded answer. The record is dropped as soon as anything else is scored, which is what keeps a refine phase, which shares the fit's trajectory, from reporting the point it started from.

When it runs

The stage is on under edition = 2 and above. Under the legacy edition it is off, because an unchanged configuration file has to keep behaving exactly as it always has, and a legacy fit with a stochastic model is told at startup which two keys turn it on.

It is also skipped in three cases. When no model is stochastic. When the wall_time_fit budget is already spent, because a budget is a promise about the whole run and this costs more simulations. And when stochastic_seed is one of the _honorbngl modes and every stochastic model pins its seed, because then every replicate would reproduce the same trajectory. Each of those says what it is skipping and why.

No fitting method searches any differently. This covers the answer the fit reports and not the search that produced it, which is the larger piece of work the issue sets aside.

Two smaller pieces this needed

make_job takes a replicate_offset. Under the default seed policy a stochastic simulation's seed comes from the parameter values and the replicate index and nothing else, so running a parameter set again at the same index reproduces the same trajectory exactly. The new stage asks for indices past every one the fit itself used. Zero, the default, is the ordinary path and is unchanged.

The scoring half of add_to_trajectory moved into score_result, so the new stage scores results it deliberately does not put in the trajectory through exactly the same path the fit used.

Verification

The full test suite passes (4788 passed, 14 skipped). 42 new tests in tests/test_best_fit_confirmation.py cover the arithmetic, the report text, the trajectory changes, the gates, and the stage itself driven by a fake dask client with a model whose objective depends on the replicate index. I checked that each of them fails when the corresponding piece of the change is removed.

For a real check, I ran the NFsim receptor example (examples/receptor_nf) for three iterations with three candidates at four replicates. The fit reported an objective of 54878. Running that same parameter set four more times gave an average of 72839 give or take 5522, so the reported value was about 18000 too good. That is the file shown above.

Not in this change

Results/information_criteria.txt computes its log-likelihood by re-simulating the best fit once, so for a stochastic model with a likelihood objective it is still a single noisy draw, and it will now disagree with the average this stage reports. Averaging a likelihood over stochastic replicates is a question of its own, so I have not touched it here. Filed separately as #676.